clear; clc; close all;
angstrom = 1e-10;  %[m]
a0 = 5.29177e-11;  %[m]
cmn1 = 0.0299792458e12; %[Hz/cm^-1]
foffset = -0.08;     %[GHz] wavemeter offset(-0.16 GHz) + AOM(+0.08GHz) 

fnames = dir('data/v6/2020*Rb2*.mat');
Nfile = length(fnames);
dataAll = cell(1, Nfile);
disp('---------------------------------------------');
disp(['There are ',num2str(Nfile), ' files in total']);
dataxAll = [];
dataRb2All = [];
dataRb2_bkg_All = [];
dataNcyc = [];

xdataAll = [];
ydataAll = [];
yerr_dataAll = [];

for i = 1:Nfile
    filename4load = [fnames(i).folder '/' fnames(i).name];
    load(filename4load);
    dataxi = xVarList1+foffset;                     %[GHz] 
    dataRb2i = xVarCounts_Rb_2;             %ion counts of Rb2
    dataRb2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
    dataNcyci = cycleNumxVar;               %cycle number
    dataxyi = [dataxi dataRb2i dataRb2i_bkg dataNcyci];
    dataAll{i} = dataxyi; 
    
    dataxAll = [dataxAll dataxi'];
    dataRb2All = [dataRb2All dataRb2i'];
    dataRb2_bkg_All = [dataRb2_bkg_All dataRb2i_bkg'];
    dataNcyc = [dataNcyc dataNcyci'];
end
fnames = dir('data/v6/more molecules/2020*Rb2*.mat');
Nfile = length(fnames);
factor3 = 2./1.63;
factor = 1.63*factor3;   %KRb number change
dataAll = cell(1, Nfile);
disp('---------------------------------------------');
disp(['There are ',num2str(Nfile), ' files in total']);
disp('---------------------------------------------');
for i = 1:Nfile
    filename4load = [fnames(i).folder '/' fnames(i).name];
    load(filename4load);
    dataxi = xVarList1+foffset;                     %[GHz] 
    dataRb2i = xVarCounts_Rb_2/factor;             %ion counts of Rb2
    dataRb2i_bkg = xVarCountsBkgd_Rb_2/factor;     %bkgnd ion counts of Rb2
    dataNcyci = cycleNumxVar;               %cycle number

    dataxyi = [dataxi dataRb2i dataRb2i_bkg dataNcyci];
    dataAll{i} = dataxyi; 
    
    dataxAll = [dataxAll dataxi'];
    dataRb2All = [dataRb2All dataRb2i'];
    dataRb2_bkg_All = [dataRb2_bkg_All dataRb2i_bkg'];
    dataNcyc = [dataNcyc dataNcyci'];
end

fnames = dir('data/v6/more molecules2/2020*Rb2*.mat');
Nfile = length(fnames);
factor = 1.42*factor3;   %KRb number change
dataAll = cell(1, Nfile);
disp('---------------------------------------------');
disp(['There are ',num2str(Nfile), ' files in total']);
disp('---------------------------------------------');
for i = 1:Nfile
    filename4load = [fnames(i).folder '/' fnames(i).name];
    load(filename4load);
    dataxi = xVarList1+foffset;                     %[GHz] 
    dataRb2i = xVarCounts_Rb_2/factor;             %ion counts of Rb2
    dataRb2i_bkg = xVarCountsBkgd_Rb_2/factor;     %bkgnd ion counts of Rb2
    dataNcyci = cycleNumxVar;               %cycle number

    dataxyi = [dataxi dataRb2i dataRb2i_bkg dataNcyci];
    dataAll{i} = dataxyi; 
    
    dataxAll = [dataxAll dataxi'];
    dataRb2All = [dataRb2All dataRb2i'];
    dataRb2_bkg_All = [dataRb2_bkg_All dataRb2i_bkg'];
    dataNcyc = [dataNcyc dataNcyci'];
end


%%%------plot-------
f0 = 447000;
fcorrection = -2.08;         %[GHz]
df = 0.2;                      %[GHz]
markersize = 5;
yplotmin = -0.05;
yplotmax = 5;

% h1 = figure('units','normalized','WindowStyle','docked');
h1 = figure('position', [0, 500, 1000, 400]);
Xplot = dataxAll- f0;
Yplot = (dataRb2All)./dataNcyc;
Yplot_err = sqrt(dataRb2All)./dataNcyc;
errorbar(Xplot, Yplot, Yplot_err, 'ok');
xlim([680,760]);
ylim([yplotmin, yplotmax]);
xlabel(['f - ', num2str(f0),' (GHz)']);
ylabel('Ion count per cycle');
title('All Rb_2 REMPI data with the rotation transitions of B^1\Pi_u(v''=6) - X^1\Sigma_g(v"=0,J")');

f_P_branch = f0 + fcorrection+[...
    687.456 692.528 697.411 702.106 706.614 710.933 715.065 719.008...
    722.763 726.330 729.710 732.901 735.904 738.719 741.346 743.785...
    746.036 748.099 749.974 751.660]; %[GHz]
f_Q_branch = f0 + fcorrection+[...
    711.034 714.982 718.743 722.316 725.700 728.897 731.906 734.726 ...
    737.359 739.803 742.060 744.128 746.008 747.701 749.205 750.521 ...
    751.650 752.590 753.342 753.906 754.282]; %[GHz]
f_R_branch = f0 + fcorrection + [755.593]; %[GHz]
col = 8;                        %plot column
NplotP = length(f_P_branch);
NplotQ = length(f_Q_branch);
NplotR = length(f_R_branch);
row = ceil(NplotP/col);          %plot row

% h2 = figure('units','normalized','WindowStyle','docked');
% for j = 1:NplotP 
%     subplot(row, col, j);
%     %%----fit---------
%     disp(' ');
%     disp(['P: j = ', num2str(22-j)]);
%     fPc = f_P_branch(j);
%     dfmax = 0.3;                      %[GHz]
%     dfmin = 0.3;
%     xfitdata = Xplot(Xplot<=fPc-f0+dfmax&Xplot>=fPc-f0-dfmin);
%     yfitdata = Yplot(Xplot<=fPc-f0+dfmax&Xplot>=fPc-f0-dfmin);
%     yfitdata_err = Yplot_err(Xplot<=fPc-f0+dfmax&Xplot>=fPc-f0-dfmin);
%     %---fit data to a model----
%     xc = 0;%fQc-f0;
%     plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
%     hold on;
%     Aguess = max(yfitdata);
%     dxguess = 0.03;
%     xfitdata = xfitdata-(fPc-f0);
%     fo = fitoptions('Method','NonlinearLeastSquares',...
%         'Lower',[0, 0, xc-0.1],...
%         'Upper',[10, 0.7, xc+0.1],...
%         'StartPoint',[Aguess, dxguess, xc]);
%     ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
%     [curvex,gofx] = fit(xfitdata', yfitdata', ft)
%     curvex1 = xc-0.6:0.001:xc+0.6;
%     plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);
%     
%     %%%-----plot-------
%     fPc = f_P_branch(j)-f0;
%     Xsubplot = Xplot-fPc;
%     errorbar(Xsubplot, Yplot, Yplot_err, 'ok', 'MarkerSize', markersize);
%     xlim([-df, +df]);
%     ylim([yplotmin, 2.4]);
%     hold on
%     plot(Xsubplot(Xsubplot<=0.17&Xsubplot>=-0.17), ...
%         Yplot(Xsubplot<=0.17&Xsubplot>=-0.17), 'ok',...
%         'MarkerSize', markersize, 'MarkerFaceColor', 'r');
%     hold off
% 
%     title(['P(J"=',num2str(22-j),'): fc = ', num2str(fPc), ' GHz'], 'FontSize', 9);
% %     xlabel(['f - ', num2str(f0),' (GHz)']);
% %     ylabel('Ion count per cycle');
% end

N_list = [];
A_list = [];
A_list_err = [];
dx_list = [];
dx_list_err = [];

% h3 = figure('units','normalized','WindowStyle','docked');
% for j = 1:NplotQ
%     row = ceil(NplotQ/col);          %plot row
%     subplot(row, col, j);
%     %%%----fit-------
%     disp(' ');
%     disp(['Q: j = ', num2str(22-j)]);
%     fQc = f_Q_branch(j);
%     dfmax = 0.3;                      %[GHz]
%     dfmin = 0.3;
%     if j==2
%         dfmax = 0.3;                      %[GHz]
%         dfmin = 0.18;
%     end
% %     if j==4
% %         dfmax = 0.18;                      %[GHz]
% %         dfmin = 0.3;
% %     end
%     xfitdata = Xplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
%     yfitdata = Yplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
%     yfitdata_err = Yplot_err(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
%     if 22-j==9
%         yfitdata = yfitdata-2.119*exp(-(xfitdata-(743.928-0.1)).^2/2/0.034^2)...
%             -3.056*(4/11)*exp(-(xfitdata-(743.928+0.027)).^2/2/0.064^2);
%     end
%     if 22-j==20
%         yfitdata = yfitdata-3.876*(14/31)*exp(-(xfitdata-(712.902+0.082)).^2/2/0.034^2);
%     end
%     if 22-j==1
%         yfitdata = yfitdata-2.06*(14/25)*exp(-(xfitdata-(752.202+0.1168)).^2/2/0.029^2);
%     end
%     %---fit data to a model----
%     xc = 0;%fQc-f0;
%     plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
%     hold on;
%     Aguess = max(yfitdata);
%     dxguess = 0.03;
%     xfitdata = xfitdata-(fQc-f0);
%     fo = fitoptions('Method','NonlinearLeastSquares',...
%         'Lower',[0, 0.04, xc-0.07],...
%         'Upper',[10, 0.07, xc+0.07],...
%         'StartPoint',[Aguess, dxguess, xc]);
% %     if 22-j == 5
% %         fo = fitoptions('Method','NonlinearLeastSquares',...
% %         'Lower',[0, 0, xc-0.07],...
% %         'Upper',[10, 0.04, xc+0.07],...
% %         'StartPoint',[Aguess, dxguess, xc]);
% %     end
%     ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
%     [curvex,gofx] = fit(xfitdata', yfitdata', ft)
% 
%     N_list = [N_list 22-j];
%     A_list = [A_list curvex.A];
%     ci = confint(curvex);
%     if 22-j==3
%         ci
%         
%         (ci(2,1)-ci(1,1))/4
%     end
%     A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
%     dx_list = [dx_list curvex.dx];
% %     if 22-j==5
% %         dx_list_err = [dx_list_err 0.001];
% %     else
%         dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
% %     end
%     
%     curvex1 = xc-0.6:0.001:xc+0.6;
%     plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);
%     
%     
%     %%%---plot------
% %     row = ceil(NplotQ/col);          %plot row
% %     subplot(row, col, j);
%     fQc = f_Q_branch(j)-f0;
%     Xsubplot = Xplot-fQc;
%     errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize',markersize,...
%         'MarkerFaceColor', 'c');
%     xlim([-df, +df]);
%     ylim([yplotmin, yplotmax]);
% %     hold on
% %     plot(Xsubplot(Xsubplot<=0.17&Xsubplot>=-0.17), ...
% %         Yplot(Xsubplot<=0.17&Xsubplot>=-0.17), 'ok',...
% %         'MarkerSize', markersize, 'MarkerFaceColor', 'c');
%     hold off
% 
%     title(['Q(J"=',num2str(22-j),'): fc = ', num2str(fQc), ' GHz'], 'FontSize', 9);
% %     xlabel(['f - ', num2str(f0),' (GHz)']);
% %     ylabel('Ion count per cycle');
% end
h3 = figure('position', [0, 500, 1000, 600]);
subplot(2, 1, 1);
row = ceil(NplotQ/col);          %plot row
for j = 2:9
    disp(' ');
    disp(['j = ', num2str(22-j)]);
    fQc = f_Q_branch(j);
    dfmax = 0.3;                      %[GHz]
    dfmin = 0.3;
    if j==3
        dfmax = 0.2;                      %[GHz]
        dfmin = 0.3;
    end
    xfitdata = Xplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata = Yplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata_err = Yplot_err(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);

    if 22-j==20
        yfitdata = yfitdata-3.876*(14/31)*exp(-(xfitdata-(712.902+0.082)).^2/2/0.034^2);
    end


    %---fit data to a model----
    xc = fQc-f0;
    plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
    hold on;
    Aguess = max(yfitdata);
    dxguess = 0.03;
    fo = fitoptions('Method','NonlinearLeastSquares',...
        'Lower',[0, 0.04, xc-0.07],...
        'Upper',[10, 0.07, xc+0.07],...
        'StartPoint',[Aguess, dxguess, xc]);
    ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
    [curvex,gofx] = fit(xfitdata', yfitdata', ft)
    
    N_list = [N_list 22-j];
    A_list = [A_list curvex.A];
    ci = confint(curvex);
    A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
    dx_list = [dx_list curvex.dx];
    dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
    
    
    curvex1 = xc-0.6:0.001:xc+0.6;
    plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);

    legend off
    hold on
    if mod(22-j,2) ==1
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', [47 108 167]/(47+108+167)*1.9, 'CapSize', 0, 'LineWidth', 1);
    else
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', 'red','CapSize', 0, 'LineWidth', 1);
    end
    
    xdataAll = [xdataAll xfitdata];
    ydataAll = [ydataAll yfitdata];
    yerr_dataAll=[yerr_dataAll yfitdata_err];
    
    ax = gca;
    ax.FontSize = 13;
    xlabel(['f - ', num2str(f0),' (GHz)']);
    ylabel('Ion count per cycle');
    xlim([712, 736]);
    xticks(712:2:736);
    ylim([yplotmin, yplotmax]);
end

subplot(2, 1, 2);
for j = 10:NplotQ
    disp(' ');
    disp(['j = ', num2str(22-j)]);
    row = ceil(NplotQ/col);          %plot row
    fQc = f_Q_branch(j);
    dfmax = 0.3;                      %[GHz]
    dfmin = 0.3;
     if 22-j==8
        dfmax = 0.1;                      %[GHz]
        dfmin = 0.3;
     end
     if 22-j==3
         dfmax = 0.3;                      %[GHz]
         dfmin = 0.1;
     end
    xfitdata = Xplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata = Yplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata_err = Yplot_err(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    
    if 22-j==9
        yfitdata = yfitdata-2.119*exp(-(xfitdata-(743.928-0.1)).^2/2/0.034^2)...
            -3.056*(4/11)*exp(-(xfitdata-(743.928+0.027)).^2/2/0.064^2);
    end
    if 22-j==1
        yfitdata = yfitdata-2.06*(14/25)*exp(-(xfitdata-(752.202+0.1168)).^2/2/0.029^2);
    end
    
    %---fit data to a model----
    xc = fQc-f0;
    plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
    hold on;
    Aguess = max(yfitdata);
    dxguess = 0.03;
    
    fo = fitoptions('Method','NonlinearLeastSquares',...
        'Lower',[0, 0.04, xc-0.1],...
        'Upper',[50, 0.07, xc+0.1],...
        'StartPoint',[Aguess, xc, dxguess]);
    ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
    [curvex,gofx] = fit(xfitdata', yfitdata', ft)

    N_list = [N_list 22-j];
    A_list = [A_list curvex.A];
    ci = confint(curvex);
    A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
    dx_list = [dx_list curvex.dx];
    dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
    
    
    curvex1 = xc-0.3:0.001:xc+0.3;
    plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);

    hold on
    legend off
    if mod(22-j,2) ==1
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', [47 108 167]/(47+108+167)*1.9, 'CapSize', 0, 'LineWidth', 1);
    else
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', 'red','CapSize', 0, 'LineWidth', 1);
    end
    
    xdataAll = [xdataAll xfitdata];
    ydataAll = [ydataAll yfitdata];
    yerr_dataAll=[yerr_dataAll yfitdata_err];
    
    ax = gca;
    ax.FontSize = 13;
    xlabel(['f - ', num2str(f0),' (GHz)']);
    ylabel('Ion count per cycle');
%     xlim([21, 40]);
    ylim([yplotmin, yplotmax]);
end

% 
% h4 = figure('units','normalized','WindowStyle','docked');
for j = 1:NplotR
% %     row = ceil(NplotR/col);          %plot row
%     subplot(row, col, j);
%     %%%----fit-------
    disp(' ');
    disp(['R: j = ', num2str(1-j)]);
    fRc = f_R_branch(j);
    dfmax = 0.3;                      %[GHz]
    dfmin = 0.3;
% %     if j==2
% %         dfmax = 0.3;                      %[GHz]
% %         dfmin = 0.18;
% %     end
% %     if j==4
% %         dfmax = 0.18;                      %[GHz]
% %         dfmin = 0.3;
% %     end
    xfitdata = Xplot(Xplot<=fRc-f0+dfmax&Xplot>=fRc-f0-dfmin);
    yfitdata = Yplot(Xplot<=fRc-f0+dfmax&Xplot>=fRc-f0-dfmin);
    yfitdata_err = Yplot_err(Xplot<=fRc-f0+dfmax&Xplot>=fRc-f0-dfmin);
    if j==1
        yfitdata = (yfitdata-5.5*(10/23)*exp(-(xfitdata-(753.513-0.053)).^2/2/0.04^2))*0.5;
    end
    %---fit data to a model----
    xc = fRc-f0;
    plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
    hold on;
    Aguess = max(yfitdata);
    dxguess = 0.03;
%     xfitdata = xfitdata-(fRc-f0);
    fo = fitoptions('Method','NonlinearLeastSquares',...
        'Lower',[0, 0.04, xc-0.03],...
        'Upper',[10, 0.07, xc+0.03],...
        'StartPoint',[Aguess, dxguess, xc]);
    ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
    [curvex,gofx] = fit(xfitdata', yfitdata', ft)

    N_list = [N_list 0];
    A_list = [A_list curvex.A];
    ci = confint(curvex);
    A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
    dx_list = [dx_list curvex.dx];
    dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
    
    curvex1 = xc-0.6:0.001:xc+0.6;
    plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);
    
    
    errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
        'MarkerFaceColor', 'r', 'CapSize', 0, 'LineWidth', 1);
    
    xdataAll = [xdataAll xfitdata];
    ydataAll = [ydataAll yfitdata];
    yerr_dataAll=[yerr_dataAll yfitdata_err];
%     %%%---plot------
% %     row = ceil(NplotQ/col);          %plot row
% %     subplot(row, col, j);
%     fRc = f_R_branch(j)-f0;
%     Xsubplot = Xplot-fRc;
%     errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize',markersize,...
%         'MarkerFaceColor', 'c');
%     xlim([-df, +df]);
%     ylim([yplotmin, yplotmax]);
% %     hold on
% %     plot(Xsubplot(Xsubplot<=0.17&Xsubplot>=-0.17), ...
% %         Yplot(Xsubplot<=0.17&Xsubplot>=-0.17), 'ok',...
% %         'MarkerSize', markersize, 'MarkerFaceColor', 'c');
%     hold off
% 
%     title(['R(J"=',num2str(22-j),'): fc = ', num2str(fRc), ' GHz'], 'FontSize', 9);
% %     xlabel(['f - ', num2str(f0),' (GHz)']);
% %     ylabel('Ion count per cycle');
end

h5 = figure('position', [0, 500, 700, 250]);
factory = 572/124;
odd = mod(N_list,2);
even = ~odd;
bar_x = N_list;
bar_y = (A_list-0.26)*sqrt(2*pi).*dx_list/0.03*factory;
bar_y_err = bar_y.*min(sqrt((A_list_err./A_list).^2),1);%+(dx_list_err./dx_list).^2
% bar(bar_x((bar_x<20)),bar_y((bar_x<20)), 'FaceColor', [0 1 1]*0.5, 'EdgeColor', [0 1 1]);
bar(bar_x((bar_x<21)),bar_y((bar_x<21)), 'FaceColor', [47 108 167]/(47+108+167)*1.3, 'EdgeColor', [47 108 167]/(47+108+167)*1.);
yticks([0:80/5:80]);
yticklabels({'0','0.2','0.4', '0.6','0.8','1'})
hold on
bar_y1=bar_y;
bar_y1((bar_x<21)&odd)=0;
bar(bar_x(bar_x<21),bar_y1(bar_x<21), 'FaceColor', [102 8 32]/(102+8+32)*0.7, 'EdgeColor', [102 8 32]/(102+8+32)*0.5);
xticks([0:20]);
xlim([-1 20.5]);
ylim([0 80]);
errorbar(bar_x(bar_x<21),bar_y(bar_x<21),bar_y_err(bar_x<21), '.', 'MarkerSize', 0.1,...
         'CapSize', 3, 'LineWidth', 1, 'MarkerEdgeColor',[1 1 1]*0,...
         'Color',[1 1 1]*0);
hold off

set(h3,'Units','Inches');
pos = get(h3,'Position');
set(h3,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h3,'plots\K2REMPI_spectrum_5G','-dpdf','-r0')

set(h5,'Units','Inches');
pos = get(h5,'Position');
set(h5,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h5,'plots\Rb2REMPI_bar_5G','-dpdf','-r0')